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We use fully self-consistent GW calculations on diamond and silicon carbide to reparametrize the 
Heyd-Scuseria-Ernzerhof exact exchange density functional for use in band structure calculations of 
semiconductors and insulators. We show that the thus modified functional is able to calculate the 
band structure of bulk Si, Ge, GaAs, and CdTe with good quantitative accuracy at a significantly 
reduced computational cost as compared to GW methods. We discuss the limitations of this func¬ 
tional in low-dimensions by calculating the band structures of single-layer hexagonal BN and M 0 S 2 , 
and by demonstrating that the diameter scaling of curvature induced band gaps in single-walled 
carbon nanotubes is still physically incorrect using our functional; we consider possible remedies to 
this problem. 


I. INTRODUCTION 

Applied materials science relies heavily on the ability 
of theoreticians to accurately predict the physical prop¬ 
erties of novel materials. Calculations are widely used as 
guidance in the design of new materials, without which 
experiments can rapidly turn into a cumbersome and ex¬ 
pensive trial-and-error process. Unfortunately, the abil¬ 
ity of theoreticians to provide guidance to experimental¬ 
ists is in practice severely limited by the accuracy of the 
methods available. One particular example for this issue 
is the problem of predicting the band gap of semiconduc¬ 
tors. 

A widely used and very accurate method in the study 
of semiconducting materials is density functional theory 
(DFT). Useful as it is, a well known problem exists with 
DFT in that it severely underestimates the band gap. 
This is especially true in low-dimensional materials such 
as boron nitride^, polyyne^*^, or even carbon nanotubes 
where even the diameter scaling of the curvature induced 
band gap is wrong, let alone its magnitude^. This prob¬ 
lem can be circumvented by going beyond the DFT level 
and using for example the Quantum Monte Carlo (QMC) 
method to calculate the band gap>^. QMC is however ex¬ 
tremely expensive to apply to semiconductors or insula¬ 
tors unless the location of the band gap in the reciprocal 
space is exactly known. 

Another approach for improving the band gap is the 
use of many-body theory, and more specifically, the GW 
approximation^. This allows one to calculate the quasi¬ 
particle band structure by starting from the DFT wave 
functions and solving the Dyson equation. The result¬ 
ing band gaps are typically in very good agreement with 
experiments and are a signifcant improvement over DFT 
results (see Figured]). The disadvantage of the GW ap¬ 
proximation is that it is computationally very expensive. 
While today it is possible to use the method on systems 
composed of as much as a hundred atoms on massively 
parallel machines^, there are plenty of calculations that 
are simply not feasible with GW to this day. Examples 


include large supercell calculations such as surface inter¬ 
calation, or the study of composite materials and het¬ 
erostructures, which are of particular interest to the sci¬ 
entific community^^— . Furthermore, in low-dimensional 
materials the GW approach converges extremely slowly 
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Figure 1: The electronic density of states for bulk diamond 
and d silicon carbide calculated with HSE06, the modified 
HSE functional, and the GW approximation. HSE06 consid¬ 
erably underestimates the gap, while our modified functional 
is in close agreement with the GW calculations. 
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as pointed out in a recent publicatioi>ii, making it all 
but infeasible for the study of nanoscale materials at this 
time. 

Due to the computational expense of both QMC and 
GW calculations, it would be advantageous to develop 
a computational method that is able to approach the 
accuracy of QMC or the GW approximation without 
the accompanying computational expenses. The triv¬ 
ial choice would be the use of exact exchange density 
functionals, also known as hybrid density functionals. 
These are known to predict a much more accurate band 
gap than traditional density functionals such as the local 
density approximation (LDA) or the semilocal Perdew- 
Burke-Ernzerhof (PBE) generalized gradient approxima¬ 
tion (GGA) as well as other GGA functionals. While 
exact exchange functionals still underestimate the band 
gap, and they in general suffer from serious problems in 
terms of predictive power as pointed out recently!^, in 
principle it is possible to modify the parameters of these 
functionals in order to deliver a better accuracy when 
calculating the band structure. 

Screened exact exchange functionals such as the Heyd- 
Scuseria-Ernzerhof (HSE) functional include dielectric 
screeningi^ which makes them ideal for use in periodic 
structures including three dimensional solids and low¬ 
dimensional materials like carbon nanotubes. These 
functionals offer two easily accessible parameters that can 
be adjusted without tampering with the core of the func¬ 
tional: these are the amount of exact exchange, and the 
screening parameter. Adjusting these parameters is pos¬ 
sible using a variety of criteria, including what has been 
dubbed the HSE12 functional^. Eor the purposes of ac¬ 
curate band structure calculations it is sensible to set 
the criteria to achieve the best possible match between 
the band structures predicted by the modified HSE func¬ 
tional and GW many-body theory. 

In this work we present a modification to the HSE func¬ 
tional with the help of fully self-consistent GW many- 
body calculations. As we show below, it is possible to ad¬ 
just the two parameters in the HSE functional such that 
the band structures calculated with the modified func¬ 
tional are in very close agreement with GW calculations. 
We compare our modified HSE functional to the HSE12 
functional. We argue that our method makes quanti¬ 
tatively accurate first-principles band structure calcula¬ 
tions on large systems feasible. The method is of course 
not without its own limitations, which we discuss at the 
end of the paper. 


II. THE MODIFIED HSE FUNCTIONAL 

In our calculations we rely on the VASP code^^, using 
a plane-wave basis set. Our goal is to develop an exact 
exchange functional that reproduces the GW quasipar¬ 
ticle band gaps in semiconductors. We use the Heyd- 
Scuseria-Ernzerhof 2006 (HSE06) functional as the basis 
for this approach, which describes the exchange part of 


the density functional as 25% exact exchange and 75% 
PBE exchange. The modified functional differs from the 
original HSE06 in two parameters only, the amount of 
exact exchange and the screening parameter. The new 
parameters are obtained by performing a least squares fit 
between the band structures obtained with our functional 
and GW band structures, then changing the parameters 
until the root mean square (RMS) is found to be minimal. 

The fitting procedure is performed using the highest 
valence and lowest conduction bands in the entire Bril- 
louin zone on a finite, 12 x 12 x 12 fc-point grid. At 
each fc-point the squared energy difference is taken be¬ 
tween the GW and the modified HSE valence band, then 
summed up over the entire grid. The process is repeated 
for the conduction band and added to the sum obtained 
from the valence band. The thus obtained RMS value 
gauges the quality of the match between the modified 
HSE functional and the GW calculation. The smaller 
the RMS the better the parameters of the modified HSE 
functional. 

The GW calculations are performed in the fully self- 
consistent approach where both the Green’s function 
and the dielectric screening are iteratively updated un¬ 
til self-consistency is reached (a maximum of four iter¬ 
ations were required). This allows for the best avail¬ 
able accuracy within GW many-body theory. Gonsider- 
ing the difficulty of converging such calculations in low¬ 
dimensional materialsii we restrict ourselves to three- 
dimensional crystals. In particular we look at diamond 
and /3 silicon carbide (3C-SiG) crystals. The geometries 
are optimized using the HSE06 functional. The GW cal¬ 
culations are then performed in fixed geometries. Both 
the structural relaxations and the GW calculations are 
carefully tested for convergence with respect to a variety 
of parameters. The convergence criteria used here are an 
uncertainty of < 0.001 A in the lattice constant during 
the HSE06 optimizations and an uncertainty of < 0.05 
eV in the energy gap during the GW calculations. The 
HSE06 optimizations are found to reach convergence at 
an 8 X 8 X 8 Monkhorst-pack grid in fc-space and 800 eV 
plane-wave cutoff energy, while the GW calculations at 
a 12 X 12 X 12 fc-point grid, 400 eV plane-wave cutoff, 16 
electronic bands, and = 48 where is the number 
of frequency points used in the contour integrations in 
the complex plane. 

With the fully converged GW band structures available 
we can make an easy comparison between many-body 
theory and existing hybrid density functionals. Eigure[T] 
shows the electronic density of states for diamond and 
/3 silicon carbide according to the GW approximation 
and the HSE06 functional. As expected the HSE06 band 
gap is much smaller. If, however, we modify the exact 
exchange ratio and screening parameters in the HSE06 
functional, it is possible to rectify this problem. 

We shall denote the two parameters as follows. Param¬ 
eter Aexx is the amount of exact exchange in the den¬ 
sity functional, i.e. the total exchange in the functional 
is Aexx ■ + (1 - Aexx) ■ where and 
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Figure 2: The root mean square of the band structure of 
bulk diamond and /3 silicon carbide between the modified HSE 
functional and GW calculations. The minimum is found at 
AEXX = 0.46 and /i = 0.24 where the RMS is 0.10 eV; 
this is a substantial improvement over HSE06 and HSE12. 


are the exact and PBE exchange terms, respec¬ 
tively. Parameter p controls the amount of screening in 
the functional; the dimension of p is inverse length and it 
determines the separation of long- and short-range terms 
in the exchange within HSE-type functionals, which is 
equivalent to applying a semiempirical screening to long- 
range interactions. In our starting point, the HSE06 
functional, Aexx = 0.25 and p = 0.2 A“^. In or¬ 
der to modify the HSE functional to fit the band gap 
to that found in GW many-body theory, we explore 
this two-dimensional parameter space within the range 
of Aexx € [0.2, 0.5] and p G [0.1,0.3] A“^. The pre¬ 
viously discussed RMS will be a function of these two 
parameters, and in order to find the ideal settings for 
Aexx and p this function must be minimized. 

Eigure[2]shows the RMS of the least square fit between 
the GW band structures and the modified HSE functional 
as a function of Aexx and p, averaged for diamond and 
/3 silicon carbide, the two benchmark materials consid¬ 
ered. The entirety of the topmost valence band and the 
bottommost conduction band is taken into account span¬ 
ning the entire Brillouin zone in a 12 x 12 x 12 grid. The 
best fit is achieved at Aexx = 0.46 and p = 0.24 A“^. 
The figure highlights the settings of the original HSE06 
functional, the HSE12 parametrization, and our modified 



AHSEmodZK 

Agw 

AhSE12 

c 

6.16 

6.05 

5.69 

Si 

1.55 

1.33 

1.35 

Ge 

1.17 

1.07 

0.89 

SiC 

2.84 

2.70 

2.50 

GaAs 

1.92 

2.12 

1.54 

GdTe 

2.14 

2.08 

1.76 


Table I: Band gaps (in eV) for the materials in Figure [3] us¬ 
ing our modified HSE functional (AHSEmodZK) and our GW 
calculations (Agvv)- Band gaps obtained using the HSE12 
functional {Ahsei 2 ) are provided for comparison. Note that 
spin-orbit coupling is neglected here, which is strictly speak¬ 
ing an important factor in GdTe. 

HSE functional. As it can be seen, the fit substantially 
improves between our modification and either HSE06 or 
HSE12. 

Now if we return to Figure [T] we can see that the den¬ 
sity of states obtained using our modified HSE functional 
with Aexx = 0.46 and p = 0.24 A“^ agrees extremely 
well with the GW calculations for both diamond and (3 
silicon carbide. 


III. APPLICATION OF THE MODIFIED HSE 
FUNCTIONAL 

Now we turn our attention to the use of the modi¬ 
fied HSE functional. First we look at its performance on 
three-dimensional semiconductors, bulk silicon (Si), ger¬ 
manium (Ge), gallium arsenide (GaAs), and cadmium 
telluride (GdTe). Figure [3] shows the band structures of 
these four materials using our modified HSE functional, 
along with those of diamond and silicon carbide. Table 
[I] compares their band gaps with our GW calculations. 
The overall agreement is quite good, with the average 
absolute deviation between the band gaps in the modi¬ 
fied HSE functional and our GW results being a modest 
0.14 eV. This is a significant improvement on the aver¬ 
age performance of the HSE12 functional which yields 
an average absolute deviation of 0.28 eV. We must note 
however that in some individual cases HSE12 performs 
better, such as bulk silicon. 

We must note here that we use geometries obtained us¬ 
ing the HSE06 functional for all our calculations. We do 
this because we cannot expect the modified functional 
to correctly predict the total energy or the Hellmann- 
Feynman forces. Density functionals are typically opti¬ 
mized to possess accurate predictive power for a broad 
range of physical properties, whereas our modified HSE 
functional is made specifically for the purposes of calcu¬ 
lating accurate band structures in a fixed geometry. This 
does not necessarily mean that the HSE06 geometries will 
be good enough, as for example the diamond lattice con¬ 
stant is notably underestimated; this is discussed below. 

Now we turn to lower dimensional materials, starting 
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Figure 3: The electronic band structure for bulk diamond (top left), bulk silicon (middle left), bulk germanium (bottom left), /3 
silicon carbide (top right), gallium arsenide (middle right), and cadmium telluride (bottom right) calculated with the modified 
HSE functional and compared to results obtained with the HSE12 functional. 
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Figure 4: The electronic band structure for single-layer hexag¬ 
onal boron nitride (top) and molybdenum disulphide (bot¬ 
tom) calculated with the modified HSE functional and com¬ 
pared to results obtained with the HSE12 functional. 



AnSEmodZK 

Agw 

AhSE12 

BN 

6.26 

7.0^ 

5.86 

M 0 S 2 

2.49 

2.84ii 

2.42 


Table II: Band gaps (in eV) for monolayer hexagonal boron 
nitride and molybdenum disulphide using our modified HSE 
functional (AHSEmodZK) and GW calculations (Agw) in the 
literature. Band gaps obtained using the HSE12 functional 
{Ahsei 2 ) are provided for comparison. 


with two-dimensional nanosheets. In particular we con¬ 
sider single layer molybdenum disulphide (M 0 S 2 ) and 
hexagonal boron nitride (h-BN). Their band structures 
are shown in Figure |4l their band gaps are compared to 
previous GW literature in Table cn 

Our band gap for M 0 S 2 is 2.49 eV which compares 
reasonably to literature estimates of 2.50 eV and 2.66 
eV quasiparticle gapsi^ obtained using the single shot 
GqWo and the partially self-consistent GWq method, re¬ 
spectively, but is notably below the fully self-consistent 
GW band gap of 2.84 eV— (albeit the experimental lower 


bound for the quasiparticle band gap is 2.5 eV— ). The 
band gap of h-BN is found to be 6.26 eV here, which 
falls between recent predictions for the band gap>i^ in 
HSE06 (5.65 eV) and single-shot GqWq calculations (7.0 
eV). While our functional significantly improves the band 
gap for both M 0 S 2 and h-BN, it still underestimates the 
GW results. This suggests that care must be taken when 
applying our functional for the study of two-dimensional 
materials. 

Finally we consider a more difficult problem, that of 
the curvature induced band gaps in carbon nanotubes. 
In single-walled carbon nanotubes (SWGNTs) the simple 
zone-folding picture predicts that whenever the difference 
between the chiral indices (n,m) is a multiple of 3, the 
nanotube is metallic. In reality, the finite curvature of the 
SWCNT surface induces a band gap in all cases except for 
armchair (n, n) SWGNTs. This means that (n, 0) zigzag 
nanotubes exhibit a small, curvature induced band gap 
when n = 3k. This has been known for many years, and 
it has been measured in STMi^i which has shown that 
the curvature induced band gap approximately scales as 
1/(P where d is the diameter of the SWGNT. Indeed, be¬ 
ing a curvature induced effect, this gap must scale with 
an even power of 1/d. DFT calculations however are 
unable to correctly reproduce this behavior. The LDA 
in fact predicts a 1/d^ scaling^, which is a clear indica¬ 
tion of the failure of semilocal density functional theory 
to accurately describe electronic excitations. Hence we 
use our modified HSE functional to recalculate the band 
structure of zigzag nanotubes in an effort to rectify this 
failing of semilocal DFT functionals. 

Figure [ 5 ] shows the band gap of zigzag nanotubes from 
(9, 0) to (20,0), with the curvature induced gaps plotted 
separately in the lower panel, using both our modified 
HSE functional and the HSE12 functional. As it can be 
seen the diameter scaling of the curvature induced band 
gap still follows the 1/d^ scaling found in LDA calcula¬ 
tions. This illustrates that the modified HSE functional 
is unfortunately not entirely suited for one-dimensional 
materials, and neither is HSE12. On the other hand, both 
functionals correctly capture the physics in the primary 
gaps of semiconducting zigzag nanotubes, as both the 
1/d scaling of the gap and the trigonal warping related 
buckling effect are correctly described. The band gaps 
are larger in our modified HSE functional, as expected. 


IV. DISCUSSION 

We have demonstrated above that empirically fitting 
the parameters of HSE-type exact exchange function¬ 
als to GW calculations is a feasible and sound approach 
to improve the performance of such functionals in band 
structure calculations. We have in particular demon¬ 
strated that by changing the parameters of the HSE06 
functional it is possible to calculate band structures that 
agree with GW calculations for the valence and conduc¬ 
tion band with an average RMS error of 0.1 eV, which is 
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Figure 5: The primary (top) and secondary (bottom) band 
gaps of SWCNTs according to our modified HSE functional 
and the HSE12 functional. The primary gaps correctly scale 
with the inverse diameter {Egap = a/d where aHSEi 2 = 9.44 
eV-A and aHSEmodZK = 9.76 eV-A; compare with the LDA 
result of aLDA = 7.1 eV-A) and the previously reported 
buckling caused by the trigonal warping effect is also ob¬ 
servable. The curvature induced secondary gaps follow a 
1/d® scaling {Egap ~ h/d^ where bHSEi 2 = 43.8 eV-A® and 
bnsEmodZK = 46.6 eV-A®; compare with the LDA result of 
bLDA = 34.1 eV-A®), which is unphysical considering that 
curvature effects should scale with an even power of 1/d. 


very impressive for a density functional. 

The importance of this finding cannot be stressed 
enough. A very simple modification enables one to use 
an HSE-type density functional to closely approach the 
accuracy of GW. Since HSE-type calculations are not 
only significantly cheaper than GW but are also much 
easier to converge, the modified HSE functional can be 
used to give accurate theoretical predictions for the band 
structure of the kind of complex materials where GW is 
simply not feasible. Therefore, a correctly parametrized 
modification of the HSE functional is invaluable for con¬ 
densed matter research, and we have demonstrated here 
that such a parametrization is possible to find. 

However, it is important to be aware of two shortcom¬ 


ings of this study. First of all, the band gaps obtained 
in our GW calculations are consistently higher than ex¬ 
pected; this is especially notable in the case of diamond 
where our GW band gap is 6.05 eV while experimentally 
it is 5.47 eV. Previous GqWo calculations performed on 
the experimental geometry yielded a much better match 
with a band gap of 5.68 eV— . Similarly, our GW band 
gaps for Si, Ge, and GaAs are notably larger than pre¬ 
vious GW calculations^!. Gonsidering that we rely on 
fully self-consistent GW rather GqWo or GWq, this is a 
surprising result. One must bear in mind, however, that 
electron-phonon coupling introduces a notable band gap 
renormalization in diamond lattices. In particular, in di¬ 
amond itself the gap is reduced by 0.4 eV— which 
accounts for most of our overestimation and leaves only 
an error of 0.18 eV unexplained. 

This remaining difference may be due to a number of 
reasons. One may be the use of HSE06 geometries in¬ 
stead of experimental lattice constants. The lattice con¬ 
stant of diamond is slightly underestimated at a value 
of 3.54 A in our calculations as opposed to the experi¬ 
mental value of 3.57 A. However, the calculated lattice 
constant of silicon is 5.44 A which matches the experi¬ 
mental 5.43 A rather well, so the lattice constant may 
not be the biggest factor. Much more likely is that our 
calculations face convergence issues. This is a surprising 
finding since, as mentioned earlier in the paper, we have 
performed a full convergence test with respect to all rele¬ 
vant parameters in the calculation, and we have done so 
without separating them (that is, they were not treated 
independently). We have used direct energy gaps at high 
symmetry points in the Brillouin zone as the convergence 
criteria. Our hnding was that these energies change less 
than 0.05 eV when we try to exceed the converged pa¬ 
rameter set. 

One possibility is that the calculations have reached 
a local minimum in the parameter space of the GW pa¬ 
rameters, and that we need to go much further for true 
convergence. In particular, it has been suggested that the 
number of empty bands should be on the order of a few 
hundred in diamond lattice materials such as silicon!!. 
On the other hand, the expectation is that the band gap 
will approach the converged value from below, while in 
our case the gap is overestimated. The reason for this dis¬ 
crepancy is currently unclear, although the next step is 
obviously to extend the convergence tests in the GW cal¬ 
culation to greater ranges in the parameter space. This, 
however, goes well beyond the scope of this work, since 
the calculations are significantly more expensive to per¬ 
form if for example the number of bands taken into ac¬ 
count is increased to a few hundred. 

We emphasize that these potential convergence prob¬ 
lems in our GW calculation do not impact the main mes¬ 
sage of our paper, namely that it is perfectly possible to 
fit an HSE-type exact exchange functional to GW calcu¬ 
lations. The parameters of Aexx = 0.46 and /r = 0.24 
A“® will no doubt need to be further adjusted before the 
modified HSE functional suggested here can be used for 
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quantitative predictions of the band gap in semiconduc¬ 
tors. In addition, in order to determine a final set of pa¬ 
rameters to use in the modified HSE functional it would 
be advisable to perform the calculations over a large set 
of materials comprising lattices of different symmetries 
and atomic composition. The parameters in this paper 
should be used only as guidelines for future calculations. 

There is a second issue with the modified HSE func¬ 
tional, that of its performance in low-dimensional mate¬ 
rials. As our GW calculations seem to overestimate the 
gap in three-dimensional semiconductors, it is somewhat 
surprising to find that the modified HSE functional un¬ 
derestimates the band gap in the case of low-dimensional 
structures. The band gaps of M 0 S 2 and monolayer boron 
nitride are both too small and the diameter scaling of 
curvature induced secondary gaps in SWCNTs is the 
same unphysical \/<P dependence as was found in LDA. 
One possible reason for this is that our functional is 
parametrized exclusively on three-dimensional crystals 
due to the difficulty in converging GW calculations in 
low dimensions. Strictly speaking the dielectric screening 
in a low-dimensional material is expected to be substan¬ 
tially different from that observed in three dimensions, 
hence in low-dimensional materials the screening param¬ 
eter /i is expected to be different. As such, the finding 
that our modified HSE functional faces difficulties in low 
dimensions is not surprising. 

The solution to the above problem is of course a 
reparametrization of the modified HSE functional for use 
in two- and one-dimensional materials. This requires a 
number of fully converged and fully self-consistent GW 
calculations to be performed on low-dimensional struc¬ 
tures to act as benchmarks. These are very difficult to 
obtain as we have seen in the case of polyyne^ and boron 
nitridei^, and our current computational resources pro¬ 
hibit us from performing them. Nevertheless the relative 
success of our modified HSE functional in three dimen¬ 
sions indicates that this reparametrization is the next log¬ 
ical step towards improved exact exchange density func¬ 
tionals for use in electronic band structure calculations. 

Finally we must note that, considering recently 
reported performance issues in exact exchange 


functional&i^, it is quite possible that such density 
functionals cannot accurately predict the curvature 
induced gaps of SWCNTs at all. The qualitative agree¬ 
ment between the HSE12 functional and our modified 
HSE functional in the case of the curvature induced 
gaps of SWCNTs further suggests that this may be 
the case. In order to convincingly answer this question 
it would be necessary to not only reparametrize the 
modified HSE functional for one-dimensional materials, 
but also to perform fully self-consistent GW calculations 
on the relevant SWCNTs, which are computationally 
highly demanding due to known convergence issues in 
low dimensions. 

V. CONCLUSIONS 

We have presented an empirical reparametrization of 
the HSE06 exact exchange density functional by fitting 
the amount of exact exchange and screening in the func¬ 
tional in order to achieve the closest match possible 
to GW calculations. We have demonstrated that this 
method is sound and accurate, as our modified HSE func¬ 
tional performs very well in three-dimensional materi¬ 
als. Since the modified HSE functional is a significantly 
cheaper method than the GW approximation, being both 
a better scaling and - more importantly - a faster con¬ 
verging method, our functional can be used to perform 
calculations on the kind of large systems that cannot be 
feasibly treated in GW while approaching the accuracy 
of GW at the same time. We have discussed the limita¬ 
tions of the modified HSE functional in low-dimensional 
materials and outlined the next steps towards the devel¬ 
opment of accurate exact exchange density functionals 
for the study of the electronic structure of nanoscale ma¬ 
terials. 
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